function model=update_sigma_02(model, MaxFun)


%global model;

N     = model.N;
k     = model.k;
k_hat = model.k_hat;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% with fmincon
%
% x0 = (model.sigma_02.^2)';
% A  = -eye(k_hat);
% B  = zeros(k_hat,1);
% 
% options   = optimset('Algorithm','interior-point','GradObj','on', 'MaxFunEvals', MaxFun, 'Display', 'off');
% 
% temp      = sqrt([fmincon(@L_sigma_02, x0, A, B, [], [], [], [], [], options, model)]');
% ind       = find(temp<model.MINVALUE);
% temp(ind) = model.MINVALUE;
% model.sigma_02 = temp; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% without fmincon

term1 = sum(model.sigma_2.^2,1);
term2 = sum((repmat(model.alpha_02,N,1)-model.alpha_2).^2,1);
term3 = ((model.eta.^2)*(sum(model.sigma_1.^2+model.alpha_1.^2,1))')';
term3 = term3*(model.lambda^2);

term4 = zeros(k_hat,1);

 for m=2:k
   for l=1:(m-1)
        term4 = term4 + 2*(model.eta(:,l).*model.eta(:,m))*sum(model.alpha_1(:,l).*model.alpha_1(:,m));
   end
 end

term4 = term4';
term4 = term4*(model.lambda^2);

term5 = 2*sum((repmat(model.alpha_02,N,1)-model.alpha_2).*(model.eta*model.alpha_1')',1);
term5 = term5*(model.lambda);

term6 = term1+term2+term3+term4+term5;

temp  = (1/N)*term6;
temp  = sqrt(temp);

ind       = find(temp<model.MINVALUE);
temp(ind) = model.MINVALUE;
model.sigma_02 = temp; 

end
